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We investigate the role of the isoscalar vector interaction and the dynamics of the Polyakov 
loop on inhomogeneous phases in the phase diagram of the two-flavor Nambu-Jona-Lasinio 
(NJL) model. Thereby we concentrate on inhomogeneous phases with a one-dimensional 
f^ ■ modulation, explicitly domain-wall solitons and, for comparison, the chiral spiral. While 

the inclusion of the Polyakov loop merely leads to quantitative changes compared to the 

^ I original NJL model, the inclusion of a repulsive vector-channel interaction has significant 

qualitative effects: Whereas for homogeneous phases the first-order phase transition gets 
weakened and eventually turns into a second-order transition or a cross-over, the domain of 
inhomogeneous phases is less affected. In particular the location of the Lifshitz point in terms 
of temperature and density is not modified. Consequently, the critical point disappears from 
^ the phase diagram and only a Lifshitz point (showing a different critical behavior) remains. 

^ ' In particular, susceptibilities remain finite. 

I. INTRODUCTION 
> 

Q^ , Due to its potentially rich and complex structure a better understanding of the phase diagram of 

, , quantum chromodynamics (QCD) still poses one of the biggest challenges in modern nuclear physics 

f — I , (for dedicated reviews see, e.g., Refs. [l|-l3|] ). Experimentally, the domain of large temperatures 

^^ . is explored by heavy-ion collisions whereas the properties at low temperatures and large densities 

are relevant for compact stellar objects. On the theoretical side ab-initio lattice calculations are 
limited to small chemical potentials and investigations at moderate densities so far mainly rely on 

^ . phenomenological models. Probably the most widely- used one in this context is the Nambu-Jona- 

Lasinio (NJL) model [J] and its extensions, sharingglobal symmetries as well as the phenomenon 
of chiral symmetry breaking with QCD (see Refs. [5|-lg] for reviews). 

In this article we explore the role of inhomogeneous phases in the phase diagram of NJL-type 
models, focusing on the inclusion of a vector-channel interaction as well as the Polyakov- loop 
dynamics. The isoscalar- vector interaction naturally arises in the NJL model when motivating 
the interaction from a one-gluon exchange in QCD SJ and was already shown to be of particular 
importance in the Walecka model at finite densities [9| . More recently, its influence on the location 
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and emergence of critical points in the phase diagram has attracted new interest [1CH12|. 
The coupling of the NJL model to the Polyakov-loop dynamics on the other hand has been in- 
troduced to mimic features of (de-)confinement, in particular at finite temperatures and vanishing 
densities |13l. Il4|. There are several open issues which complicate the setup of the model at finite 
chemical potential [15|, [la]- However, at vanishing temperatures, the Polyakov loop always decou- 
ples from the quark sector by construction. In the regime of low temperatures and (not too) high 



chemical potential one typically finds a chirally restored but "confined" phase, which is sometimes 
related to the quarkyonic phase [iTj], suggested in Ref. ISj] for QCD in the limit of a large number 
of colors (Nc). 

Although most studies are usually restricted to homogeneous phases, the importance of chiral 
crystalline phases being characterized by an inhomogeneous order parameter has been pointed out 
long ago. Well known examples are the Skyrme crystal 19(] and the chiral density wave (CDW) 20]. 
For quark matter, chiral crystalline phases have been explored in the weakly coupled regime for 



large Nc J2ll. |22|]. where they form the ground state at vanishing temperatures even at asymptotic 



densities. Related to that they naturally show up in the quarkyonic matter picture 23] and in 
holographic models [2j]. For metals, however, the underlying mechanism of particle-hole pairing 
was however already considered in the 1960s 25|]. In a similar context inhomogeneous phases have 
also been discussed for color sup er conductivity, see e.g. Refs. 2614321] . again borrowing ideas from 
condensed matter physics [33|, |34 ] . 

An inhomogeneous phase also exists in the large A'^c-limit of the 1 + 1-dimensional Gross-Neveu 
(GN) model 35(, where one has a rather complete picture of the phase diagram |36l439l] . Here the 
chiral crystalline phase emerges at low densities by formation of kink and antikink solitons and 
reaches out to arbitrarily high densities for sufficiently small temperatures. 

In the NJL model, most investigations of inhomogeneous phases have been performed for the CDW 
("chiral spiral"), corresponding to a single plane wave ansatz for the order parameter [40|, |41|. The 
notable exception here is Ref. 4^. More recently, it was shown that the known solutions from 
the 1 + 1-dimensional GN model can be used to construct solutions of the 3-1-1 dimensional NJL 
model whose order parameter varies in one spatial direction 43|] . These solutions are more favored 
than the chiral spirals and the corresponding phase occupies a larger region of the phase diagram. 
However, due to the short-ranged interaction (in contrast to the weak-coupling QCD analysis at 
large Nc) and different kinematics compared to the GN model, this phase is constrained to a finite 
range of densities. 

In the chiral limit, we can thus distinguish three different phases, namely the homogeneous chirally 
broken phase, the inhomogeneous phase, and the chirally restored phase. These phases meet in a 
single point, a so-called "Lifshitz point" (LP). Within a Ginzburg-Landau analysis it can be shown 
that in the NJL model without vector interactions the LP exactly coincides with the (tri-) critical 
point (CP) of the chiral phase transition, which would be present if the analysis was limited homo- 
geneous phases 4J]. In Ref. 43| this was confirmed by an explicit model calculation. Moreover, it 
was found that the inhomogeneous phase is bordered by second-order phase boundaries, and that 
it completely covers the would-be first-order transition line between the homogeneous phases. 
In the present work we extend these investigations to include repulsive vector interactions. As 
the main result we find that the LP stays at the same temperature, whereas, when the analysis is 
restricted to homogeneous phases, the CP is shifted to smaller temperatures and larger chemical 
potentials. The critical region surrounding the CP thus disappears from the phase diagram as 
it is replaced by an energetically more preferred inhomogeneous ground state. Consequently, the 
divergence of susceptibilities at the CP, which has been related to event-by-event fluctuations in 
heavy-ion collisions and discussed as one of the most important experimental signatures of the 



CP [45|, is removed from the phase diagram. For completeness, we also discuss the inclusion of the 
Polyakov-loop dynamics, which affects the structure of the phase diagram only quantitatively. 
This paper is organized as follows: In section |TI] we discuss the role of the isoscalar-vector in- 
teraction when investigating the NJL model's phase diagram allowing for inhomogeneous phases. 
After introducing the model and the general framework for investigating inhomogeneous phases in 
subsections III Al and IIIBl respectively, we move on to an extensive numerical study of the phase 
diagram in subsection III CI The latter includes results for the density profiles in solitonic ground 
states, a comparison to the CDW, a discussion of number susceptibilities and the role of finite 
current quark masses. Part of the obtained results is then explained in the context of a generalized 
Ginzburg-Landau (GL) expansion in subsection III Dl The subsequent section IIIIj is devoted to the 
inclusion of Polyakov loop dynamics and its effect on the structure of the phase diagram. Finally, 
we summarize and conclude in section HVl 

II. THE TWO-FLAVOR NJL MODEL WITH VECTOR-CHANNEL INTERACTION 

A. Model and Approximations 

In view of finite density investigations, an important extension of the original NJL model is given 
by the inclusion of an isoscalar vector-channel interaction. The Lagrangian is then given by 

C = ij {i^^'^^ - m) V + G5 ( (^V) ' + (V^^tV V) ') - Gv {^J^i^) ' , (1) 

where ^ is a 4A'^jA'^c-dimensional quark spinor for A^^ = 2 flavors and Nc = 3 colors, j^ and t"" are 
Dirac and Pauli matrices, respectively, and m is the degenerate current quark mass. 
In mean-field approximation we expand around the expectation values of the non-vanishing, 
potentially spatially varying scalar condensate S'(x) = (^(x)^(x)), pseudo-scalar condensate 
P3(x) = ('0(x)i7^T^i/;(x)) and density n(x) = (^'''(x)^(x)), yielding 

{iPi^f ~ -S(x)2 + 25(x)V5^ , 
{iPij^T^Tpf ~ -P(x)2 + 2P(x)^i7V3V' , 

{iljjf'iljf ~ -n(x)2 + 2n(x)^7°V- (2) 

In the case of a periodic condensate with Wigner-Seitz cell V and using the imaginary-time for- 
malism, the mean-field thermodynamic potential is then given by 

n = -77 In / VlpVll; ex.p / (£mean-field + /^■!/'7°'0) 

^ kinetic ~r ^ 'cond j \y) 
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the quasi-particle Hamiltonian 

H-n = -i-i^fdi + \- (M(x) + M(x)* + 7V3m(x) - 7V3m(x)*) -/i(x) , (5) 



the spatially dependent "constituent quark" mass 

M(x) = m-2Gs(5(x)+iP(x)), (6) 

the renormalized quark chemical potential 

/i(x)=/x-2Gyn(x), (7) 

and Matsubara frequencies u^ = (2n + 1)ttT. 

For general spatially dependent mean fields, the evaluation of the functional trace in Eq. ^ and 
the subsequent minimization of Vt is highly non-trivial. However, as shown in Ref. [43], for the 
case without vector interactions a solution can be found if only one-dimensional mass modulations 
are considered. In the present work we want to generalize these solutions to Gy > 0. Since this 
cannot be done exactly in a straightforward manner, we approximate the density in Eq. ([7]) by its 
spatial average, 

n(x) — )■ n = (n(x)) = const. (8) 

As a consequence, /i becomes constant as well, and the problem reduces to the known case without 
vector interaction at a shifted value of the chemical potential. Of course, at first sight, it seems 
rather questionable whether the replacement Eq. (J16p is a good approximation in an inhomogeneous 
phase. However, as we will discuss in more detail later on, it can be rigorously justified in the 
vicinity of a second-order phase boundary to the restored phase, and in particular for the Lifshitz 
point. 

The kinetic contribution Okjnetic to the thermodynamic potential can now be evaluated from the 
eigenvalue spectrum {Ei} of the Hamiltonian Hq through the relation 

T^lnl^^iicJn + E,)^ = ^E, + T\n(^l + exp(^-^^^ . (9) 

Introducing a density of states to express — w X]_e ~^ —'^^cfdEp^E), we then formally obtain 
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-2N,jdEp{E)f{E), (10) 



with 



f{E) = /uv(-S) + /thermal (-E) , 

fvv{E) = E, 
/thermai(i^) = T In (l + exp (- ^^\\ + T In (l + exp (- ^±^) ] , (11) 



where we have used the fact that the eigenvalues of Hq come in pairs {Ei, —Ei}. 
The last missing step is a regularization of the diverging integration. For homogeneous phases 
this is mostly done by a momentum cutoff (see, e.g., Refs. [5|48|]). This is however not possible for 
inhomogeneous phases, since the quasi-particle energies can no longer be labelled by a conserved 
three-momentum. Instead we have to apply a regularization of the functional logarithm, e.g., 
by a proper-time regularization, which essentially acts on the energy spectrum rather than the 
quasi-particle momenta. With the density of states growing like p{E) = E'^/vr^ -|- 0{E^) in three 
dimensions, we then only have to regularize the contribution at T = /I = stemming from /uv (-£')• 



We regularize this part by a specific b^ 
Pauli-Villars regularization of the form g] 



ocking function in the proper-time integral leading to a 



fvv{E) ^ fME) = Y,c,^/E'^+jk\ (12) 

i=o 

with Co = 1, ci = —3, C2 = 3, C3 = —1 and a cutoff scale A. 

Within this setup we aim to discuss ground states by considering the stationary conditions 

5Vt 



dM{x) 



0, (13) 



determining the order parameter M(x) and the shifted chemical potential jl. We recall that the 
latter is a space- independent constant after inserting our approximation Eq. ([8]) into Eq. (j7]). For 
the later discussion it is helpful to cast condition p^ into 

fi = /2 + 2Gyn, (15) 

where the average density is given by 



n = 2Nc I dE p{E){n+ - n^) (16) 



with the occupation numbers n± = 1/(1 + exp((i5^=F/z)/T)). We emphasize that, besides explicitly 
depending on /2 via the occupation numbers, n is a functional of M(x) via the density of states 
p{E). 

At this stage it is worth noting that for most quantities all dependence on Gy can be absorbed 
into fl. In particular the form of the gap equation (J13p and, hence, its solutions M(x) are identical 
to the case without vector interaction upon replacing /x — )• /i. As obvious from Eq. ()16p . the same 
is true for the average density n. As a consequence, the mass functions M(x) at a given n do not 
depend on Gy- For homogeneous phases this is a well-known result [g]. The remaining effect of Gy 
is on the one hand side to map p onto /x via Eq. (J15p . and on the other hand to shift the value of 
the thermodynamic potential of a solution by —Gyf)?. This will be important for the explanation 
of our results later on. 



B. One-dimensional modulations 



For the explicit evaluation of the expressions derived in Sec. IIIAI we still need to determine the 
density of states p{E). In general, for arbitrary mass functions M(x), this is highly non-trivial 
and the problem gets even more involved for the subsequent variation in M(x). However, as 
mentioned earlier, if we restrict ourselves to phases with one-dimensional modulations the task can 
be reduced to a problem in the 1 -|- 1-dimensional GN model [4^. Since furthermore for the GN 
model all inhomogeneous phases have been classified [4f|] and in particular the phase diagram has 
been discussed in detail [36|439|, the investigation of these phases simplifies strongly. 
In the following we always assume that the one-dimensional modulations are in z-direction, i.e., 
the system is invariant under translations in the xy-directions. In the chiral limit the favored mass 
functions then take the form^ 

sn(Az|z^)cn(Az|z/) 



M,oliton(x) = z^A 



(17) 



dn(A2;|z^) 

where sn, en, and dn are Jacobi elliptic functions, and u and A are two independent parameters. 
-^soiiton(x) parameterizes a lattice of domain- wall solitons. For zv — )• 1 this solution becomes 
thermodynamically degenerate with a homogeneous phase, since this limit corresponds to a single 
soliton localized around z = 0, which does not contribute in the thermodynamic limit. 
The density of states in this phase is explicitly given by 43 1 
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where K and F are the complete and incomplete elliptic integrals of 1st kind, respectively, and 

E are the (complete or incomplete) elliptic integrals of 2nd kind. Furthermore we introduced the 

notations u = 1 — u, 9 = arcsin(£'/(\/^A)), and 6 = arcsin(A/£'). 

With these expressions at hand, the stationary condition (J13p for the ground state of the system 

is reduced from a complicated functional derivative to the much simpler problem of extremizing 

the thermodynamic potential in the two parameters A and v. As argued in Ref. 43| the solutions 

obtained by extremizing the thermodynamic potential in the remaining parameters then also fulfill 

the more general stationary condition (|13p . 

We can also consider finite current quark masses, for which the order parameter is generalized to 

cn(6|z^)dn(6|i/)~ 



Mioiiton,m(2) = A i/sn(6|i/)sn(A2:|z^)sn(Az + 6|i/) + 



sn(6|z^) 



(19) 



^ This expression can be rewritten as Msoiiton(x) = \^A' sn{A' z\u') where u' = (i- 
but Eq. lfT7)l is more convenient for our purpose. 



'and A' = (l+%/r~^)A, 



Here b is an additional parameter to be varied together with A and z/. The density of states changes 
to 

E 



where 6 G [0, oo] is given by sn(6|z^) = ,^ . The chiral limit, i.e. m = 0, corresponds to (5 = or 
equivalently b = K.{i>). 

For comparison we will sometimes also consider the chiral density wave (CDW or "chiral spiral") 
defined by the mass function 

McDw(x) = A exp{iqz) . (21) 

In this case we restrict ourselves to the chiral limit and the corresponding density of states is given 



by|43| 
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PCDwiE) = ^{ 9{E-q- A)V(i?-g)2-A2 



+ 9{E-q + A)e{E + q - A) y^ {E + q)^ - A^ 

-A-E) ( V(i? + q? - A2 - y'iE - qf " A2) } . (22) 



In the NJL model, however, the CDW is always less favored than the soliton lattice. The main rea- 
son to include it in our discussion is its simplicity. In particular its density profile is uniform. This 
is most easily seen by applying a global chiral transformation of the form ip — )• exp(i75T3 qzQ/2)^ 
with some constant zq. While the density {'ip'{'K)'tp{'x.)) is invariant under this transformation, it 
leads to a phase shift in M(x) which is equivalent to a translation by zq in the z-direction. Hence 
the density must be uniform. The replacement Eq. ^ is therefore an exact manipulation for the 
CDW, which will allow us to comment on the corresponding approximation for the solitons. 

C. Numerical results including the vector-channel interaction 

Having set up the formalism for the inclusion of the vector interaction, we now turn to the numerical 
part of our investigations. As in Ref. [43[] we fix the parameters in the chiral limit by requiring for 
the pion decay constant /tt = 88 MeV and for the constituent mass in the vacuum Mq = 300 MeV. 
The resulting parameter values are A = 757.0 MeV and Gs = 6.002/A2. The value of Gy, on the 
other hand, is treated as a free parameter, which will be varied in order to study its effect on the 
phase diagram. Starting from a color-current interaction and performing a Fierz transformation 
one would get Gy /Gs = 1/2, whereas fits to the vector meson spectrum typically yield larger 
values [a, |47|] . In the following we will therefore vary Gy between zero and Gg ■ In particular we 
will restrict ourselves to repulsive vector interactions, similar as in the Walecka model. 

1. Phase diagrams in the chiral limit 

Probably the most interesting result of this work is the effect of the vector-channel interaction 
on the transition lines in the phase diagram. In Fig. [1] we present the fJ- — T phase diagrams for 
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FIG. 1; The phase diagram in the chiral limit for different values of the vector coupling when allowing for 
the domain-wall soliton lattice. The black dashed lines represent the second-order transition lines joining 
at the Lifshitz point (dot), the shaded region represents the inhomogeneous phase. The blue solid lines 
represent the first-order phase transition obtained when limiting to homogeneous order parameters, which 
turns to second order (blue dot-dashed lines) at the critical point (square). 



various values of Gy, focusing on the region where the domain- wall soliton lattice is preferred. 
The calculations have been performed in the chiral limit. For comparison we have also indicated 
the transition lines one obtains when the analysis is hmited to homogeneous phases. In this case 
there is a critical point for small values of Gy (blue square) below which the transition from the 
broken to the restored phase is first order (blue solid line). Upon increasing Gy the CP is shifted 
to smaller temperatures (and higher chemical potentials) and eventually hits the zero temperature 
axis, so that the first-order phase transition is absent for larger values of Gy- This behavior is 
well known 48|, |49|] and has recently attracted new interest in the discussion of the critical surface 

However, the picture changes considerably, when we allow for inhomogeneous solutions. We then 
always find a regime where the domain- wall solitons are preferred (shaded region), so that we 
can distinguish three different phases: the homogeneous broken phase, the restored phase and 
the inhomogeneous phase. The corresponding three phase boundaries are all of second order. 
Their conjunction defines a Lifshitz point (dot), above which the phase boundary coincides with 
that found when limiting to homogeneous phases. Moreover, for Gy = the LP precisely agrees 



with the CP of the purely homogeneous analysis [43|, |4j]. It turns out, however, that this is 
no longer true for Gy > 0: Whereas with increasing vector coupling the CP moves downwards 
in temperature and eventually disappears from the phase diagram, we observe that the Lifshitz 
point is only shifted in the ;U-direction, while remaining at the same temperature. Consequently, 
unlike the first-order boundary in the purely homogeneous case, the existence of the inhomogeneous 
phase is not inhibited by the vector interaction. At vanishing temperature the transition from the 
homogeneous broken to the inhomogeneous phase is only slightly varying with Gy, whereas the 
transition from the inhomogeneous to the restored phase significantly shifts, thus enhancing the 
domain where inhomogeneous phases are favored. 

The observed Gy-dependence of the phase diagram can easily be understood when we recall from 
the end of section III Al that the stationary condition (J13p depends on Gy only through /i. Conse- 
quently its solutions M(x) at given (/i, T) are independent of Gy and thus equal to the solutions 
at Gy = where jl = fi. At Gy > 0, we can then translate these solutions into solutions at 
shifted values of /i, which are obtained from Eq. (jlSp . For a given mass function this mapping is 
unique since the average density n, which enters Eq. (fT5]) . also depends on Gy only through /2, see 
Eq. dm). 

The consequences of this mapping can be elaborated further for the phase transition lines. We first 
consider the case of second-order phase transitions as found in our calculation. When approaching 
the second-order transition from any direction, for each extremum in the thermodynamic potential 
that is relevant for the transition the density approaches the same value and consequently all 
extrema are mapped to the same value of //. Since there is only one extremum left on one side of 
the phase transition, we only have to make sure that the additional extremum on the other side 
cannot be mapped beyond the phase transition line, thus generating a spinodal region. For the 
transitions into the restored phase this is guaranteed by the fact that the restored solution has 
the highest possible density at given fl and that the density increases with jl. Therefore no other 
solution can be mapped beyond the transition line for Gy > 0. Similarly the homogeneous broken 
solution has the lowest density at given jl and therefore no inhomogeneous solution can be mapped 
below the transition line when going from the inhomogeneous into the homogeneous broken phase. 
Consequently all second-order transition lines for Gy = are mapped onto second-order transition 
lines for Gy > 0. In particular, this explains why the Lifshitz point stays at the same temperature, 
as the mapping leaves T untouched. 

From the arguments given above, it follows immediately that the phase diagram in the fl-T plane is 
independent of Gy. Moreover, since in this case n{T, jl) is uniquely given by Eq. (|16p and therefore 
also independent of Gy, this means that the n-T phase diagram is independent of Gy as well. This 
diagram is presented in Fig. [2j 

A slight complication of this picture arises if there is a first-order phase transition at Gy = 0, which 
occurs when limiting to homogeneous phases. In this case we have a spinodal region enclosing the 
first-order phase transition line, where the thermodynamic potential has several extrema at the 
same value of //: two local minima and one local maximum. As before, this means that at Gy > 
we have several solutions at the same value of jl. However, since the local extrema correspond to 
different masses and therefore to different densities, they will now be mapped onto different values 
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FIG. 2: The phase diagram m the n ~~ T plane. The transition hnes do not depend on Gv for Gv > 0. 



of fi by Eq. (jlSp : More precisely, solutions with lower masses will be shifted to higher values of ^ 
than solutions with higher masses. As a consequence, the spinodal region shrinks with increasing 
Gv, i.e., at fixed temperature the first-order transition gets weakened and eventually becomes 
second order. This explains why the CP moves to lower temperatures and finally disappears from 
the phase diagram. (See also Ref. [11(] for a detailed discussion of this effect.) 
Finally, we would like to comment on the slopes of the three phase transition lines in the LP. For 
Gv = all of them are equal, i.e., the phase boundaries are tangential in the LP. In the /2 — T 
plane, this remains of course true for Gv > 0. However, when we employ Eq. (|15p to map fl onto 
II we have to keep in mind that the average density n at given T and jl depends on the constituent 
quark mass M. The latter vanishes identically at the two phase boundaries to the restored phase, 
but is nonzero at the boundary between the homogeneous broken and the inhomogeneous phase, 
approaching zero only towards the LP. Depending on the corresponding critical exponent, this can 
affect the slope of this boundary in such a way that it meets the two others with a nonvanishing 
angle at Gv > 0. As one can see in Fig. [H this is indeed the case. 



2. Density profiles 



The results presented above are based on the assumption that in the thermodynamic potential the 
density n(x) can be approximated by its spatial average n, see Eq. ([8]). Of course this approximation 
might appear questionable in an inhomogeneous phase. In order to get some insight into this issue, 
we now want to discuss the density profile n{z) of the solitons. To that end we restrict ourselves 
to the case Gv = 0, where Eq. ^ does not enter into the derivation. To the extent that Eq. ([8]) 
is a good approximation, the results can then be translated to Gv > by the mapping Eq. (jlSp . 
For simplicity, we will also limit ourselves to the chiral limit, although this is not essential. 
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In mean-field the density profile is given by the expectation value 



n(x) = (V^(x)V^(x)) 

= ^ E^k W^i^. W (^+(^*) - ^-(^^)) ' (23) 



2^.. 



where the ■0Ei(x) are eigenfunctions of the Haniiltonian Hq (see Eq. ([5])) for the eigenvalues Ei and 

n±{E) are the Fermi occupation numbers defined below Eq. (J16p . 

For one-dimensional modulations the energy spectrum is highly degenerate and we can label the 



eigenvalues as E = a/A^ + Pj^, where p_L is the conserved momentum of the quasi-particle per- 
pendicular to the modulation and A is the eigenvalue of the Hamiltonian at p_L = 0. The latter 
reduces to the Gross-Neveau Hamiltonian and we refer to Ref. 43| for details. Since ip^{z)ilj^{z) 
in this case only depends on A, we arrive at 

n{z) = ^ f dXp,{\) f -^4{z)i;,{z){n+{E)-n^{E)) , (24) 



where pi{X) is the spectral density of the one-dimensional GN model 36l439l| and 



-,t,.w, u^ - (A/A)^ + i((M(z)/A)2 + . - 2) 



^^(^)^^(^) = (a/aV-e(.)/k(.) • (2^) 



Similar to the determination of the effective density of states, Eq. ([TH]) (see Ref. [43] for details), 
it is then a tedious, but straightforward exercise to cast the expression for the density profile into 
the form 

nsoiiton(^) = 2N, / dEpD,sohtoniE,z){n+{E)-n^{E)) , (26) 

where the density matrix element PD,soiitoniE, z) can be related to Psoiiton(-£')7 Eq. (fTSl) . upon the 
replacement 

PD,sohton{E,z) = Psoliton(-E) ^(^ „ i /"/ Mil^ ^ ■ "^' ^'^'^^ 



The resulting density profiles at T = and four different chemical potentials are shown in the lower 
part of Fig. [3l Comparing them with the corresponding mass functions M(z), which are displayed 
in the upper part of the figure, we find that the density distributions follow closely the positions of 
the solitons, i.e., of the zero-crossings of the mass functions. In a bag-model like picture, this can 
be interpreted as the quarks being squeezed by the bag pressure of the domains with broken chiral 
symmetry into the regions of space where chiral symmetry is almost restored. From a topological 
point of view, these quarks are related to zero energy modes localized on the domain-wall. 
These features are seen most clearly at /i = 307.5 MeV, which is just above the phase boundary 
from the homogeneous broken phase. Here the solitons are well separated, leading to strongly 
localized density peaks as functions of z. In this regime the assuption of a homogeneous density 
in order to obtain the solutions for Gy > is certainly poor. However, when p is increased the 
solitons quickly start to overlap and the density profiles become more and more washed out. At 
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FIG. 3: Spatially dependent constituent mass function M(z) (upper row) and corresponding density profile 
(lower row) at Gy = and m = for T = and, from left to right, fi = 307.5, 309, 325, and 345 MeV. The 
density is normalized to the density in the restored phase, Urest — 3— f M'^- 



the second-order transition to the restored phase, the order parameter melts and the density profile 
smoothly approaches the uniform density of the restored phase. This is illustrated by the example 
oi jJL = 345 MeV, which is close to the phase boundary. Note, however, that already at ;U = 325 MeV 
the density variations are relatively small. 

The T = results shown in Fig. [3] are the most extreme cases. With increasing temperature the 
mass functions melt and the density profiles become washed out by thermal effects. The assumption 
of a uniform density profile is therefore most questionable at the transition from the homogeneous 
broken to the inhomogeneous phase at low temperatures. In fact, at Gy > 0, a local density 
approximation would result in a lower value of fi within the solitons. This reflects the repulsive 
nature of the vector interaction, which disfavors localized density peaks. We therefore expect 
that the formation of well separated solitons, as we find near the boundary to the homogeneous 
broken phase at T = 0, eventually becomes inhibited by the vector interaction. Related to this, 
the second-order phase transition from the homogeneous broken to the inhomogeneous phase may 
partially turn into a first-order transition at Gy > 0. 

On the other hand our assumption of a uniform density becomes gradually better with increasing 
/i or T and is fully justified at the phase transition line to the restored phase. In particular the 
phase boundary itself and the Lifshitz point are not affected by the approximation. This will be 
shown more rigorously in Sec. IIIDl bv a Ginzburg-Landau analysis. 

3. C'hiral density wave 



Complementary support for the results of Sec. IIIC II can be obtained from investigations of the 
chiral spiral, Eq. (|2ip . Although, at least at Gy = 0, this kind of modulation is disfavored compared 



to the domain- wall soliton 43|], a constant density profile is no assumption here, but a property 
of the state. We can therefore perform a completely self-consistent mean-field calculation, leading 
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FIG. 4: Phase diagrams in the fJ. ~ T plane at Gv — (left) and Gv — Gs (right). The dashed lines have 
the same meaning as in Fig. [T] and represent the second-order phase boundaries between the homogeneous 
broken phase, the restored phase and the inhomogeneous phase with domain-wall solitons, respectively. The 
shaded area indicates the region where the CDW is favored compared to homogeneous phases (but not to the 
solitons). Here the blue solid line represents the first-order phase boundary from the broken homogeneous 
phase to the CDW, while the phase boundary from the CDW to the restored phase is second order and 
coincides with the boundary from the soliton phase to the restored phase. 



to the phase diagrams shown in Fig. [H The calculations have again been performed in the chiral 
limit and for Gy = (left panel) and Gy = Gs (right panel). The regions where the CDW is 
favored against the homogeneous broken and restored phases are indicated by the shaded areas. 
For comparison we have also indicated the boundaries of the regime where the solitons are favored 
(dashed lines). 

As discussed in section IIIDI below, the transition from the chiral spiral to the chirally restored 
phase is second order and agrees exactly with the transition from the soliton lattice to the restored 
phase. In particular, this also holds for the Lifshitz point. On the other hand, the transition from 
the homogeneous broken phase to the state where the chiral spiral is preferred is first order. For 
this reason, given the arguments of Sec. IIIC 1[ it is directly depending on Gy, not only through /I. 
As a result, the phase boundary, which at Gy = almost coincides with the corresponding phase 
boundary of the soliton phase, moves away from it at Gy = Gs- However, this effect is rather 
mild. Moreover we find that the first-order transition line from the homogeneous broken phase to 
the chiral spiral is much less affected by the vector interaction than the second-order transition 
line from the chiral spiral to the restored phase. The qualitative behavior of the phase diagram 
as a function of Gy is therefore similar to our results for the soliton lattice. Since we expect the 
CDW to be disfavored compared to the solitons, we can take the CDW result as a lower limit for 
the area occupied by an inhomogeneous phase. This suggests that the effect of assuming a uniform 
density profile for the solitons is not drastic. 
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4- susceptibilites 

The divergence of susceptibilities near a critical point has led to suggestions for how to locate the 
CP in the QCD phase diagram experimentally [45( and therefore attracted significant interest also 



in model studies. Since the CP as a cornerstone of the phase diagram is essentially replaced by the 
LP in our study, we first want to discuss the behavior in its vicinity. For simplicity we will limit 
ourselves to the number susceptibility 

d'^Q dn . , 

Xnn = -^ = ^, (28) 

which corresponds to the change in density when going along a line of constant temperature. For 
this reason the behavior of Xnn near the LP is in fact determined by the behavior when going from 
the homogeneous broken to the restored phase and not by the inhomogeneous phase. Consequently, 
since the LP coincides with the CP for Gy = 0, we find the same divergent behavior as when limiting 
to homogenous phases for Gy = 0. For completeness, this is illustrated on the right hand side of 
Fig. [5l showing a l/-y//icr — /i-like singulartity when approaching the LP from the left. However, 
for Gy > the behavior is qualitatively altered: The would-be CP when limiting to homogeneous 
phases is hidden inside the domain of a inhomogeneous phase and the LP does not correspond to 
the endpoint of the second-order phase transition line when limiting to homogeneous phases. For 
this reason the number susceptibility does not diverge near the LP for Gy > 0. This can easily be 
understood in the context of a Ginzburg-Landau expansion as introduced in subsection III Dl and is 
related to the fact that C4^a 7^ for Gy > at the LP. 

Focusing on Gy = first and investigating the number susceptibility numerically in the regime 
where soliton lattices are energetically preferred, we find the results shown in Figs. O [6l Since 
all transition are second order, the averaged number density changes continuously when varying 
temperature and chemical potential. At the transition from the homogeneous broken to the inho- 
mogeneous phase the change is however very rapid as discussed in the following and Xnn diverges. 
This is most prominent at T = and decreases when going towards the LP. 

In order to get an intuition for the qualitative behavior of Xnn we consider the ON model and focus 
on the density as a function of chemical potential nGjv(/i) at vanishing temperatures. As can be 
extracted from Refs. 



36l . l38l | it is given by 



ir^ 



where Mq as the fermion mass in the vacuum sets the scale and the elliptic modulus as a function 
of chemical potential is given through the implicit relation vry^/i = 2E(i/)Afo- At the transition 
from homogeneous to inhomogeneous phase at ficr = -Mq the number density then behaves like 

ln{fl/flcr - 1) ' 

which leads to a {{fj, — ficr) log^(/U — /icr)]~^-like singularity in the number susceptibility. 
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FIG. 5: Averaged density (upper panel) and number susceptibility (lower panel) for vanishing temperatures 
(left) and for the temperature at the LP (right) as a function of chemical potential and for Gv — 0. 



A straightforward generalization of Eq. IpO]) to three spatial dimensions suggests that for Gy = 
the change An in the average number density near the transition from the broken homogeneous to 
the inhomogeneous phase at /i = Hcr(T) should be given by 



An 



cn^j. 



ln{fi/fic 



1 



(31) 



with some temperature-dependent coefficient c. Indeed, our numerical results are consistent with 
this behavior, thus explaining the divergence of Xnn at ^ = fj-cr- 

In contrast, for Gy > the mapping /i — 5- /u via Eq. lfTS]) leads to a qualitatively different behavior. 
For this case we find 



Xn 



dn d^ 
djl djl 



1 



dn 



l + 2GyfS5/l 



(32) 



A(a«) 



Therefore a divergence in dn/d^ at Gy = does not result in a divergent number susceptibility 
for Gy > 0, but merely leads to a jump of order 1/2^^. This is illustrated on the right hand side 
of Fig. El 
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FIG. 7: /i — T phase diagram for different values of Gy and a current quark mass m = 5 MeV. The dots 
indicate the Lifshitz points above which the second order transition lines join and turn into a crossover. 



5. Finite current quark masses 

For completeness we present in Fig. [7] our results for the phase diagram at a non- vanishing current 
quark mass m = 5 MeV and various values of Gy by considering Eq. ()19p . Since there is no exact 
order parameter to distinguish the homogeneous broken from the restored phase in this case, there 
are strictly speaking only two phases ~ a homogeneous and an inhomogeneous phase - and, hence, 
no Lifshitz point. Nevertheless we can easily identify the remnant of the LP as a cusp in the phase 
boundary. For simplicity, we will call this point a Lifshitz point as well. 

It has already been observed |43| that the LP shifts to smaller temperatures and larger chemical 
potentials when increasing m. The dependence on Gy however stays qualitatively the same as 
in the chiral limit: the LP is only shifted in the //-direction upon increasing Gy and the domain 
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of inhomogeneous phases in the ^J- — T diagram is stretched. The explanation for this behavior is 
identical to the one given in section [lIC II for the chiral limit. 

D. Ginzburg-Landau expansion 

In the vicinity of a second-order phase transition and in particular of a critical point, the thermody- 
namic potential can be studied systematically within a Ginzburg-Landau expansion. In the present 
case, this corresponds to an expansion of the thermodynamic potential as an effective action in 
(5M(x) = M(x) — Mq and 5/u(x) = /i(x) — /io around their values M(x) = Mq and 6p,{x) = p,Q in 
the restored phase. For simplicity, we restrict ourselves to the chiral limit, so that Mq = and 
therefore 5M(x) = M(x). However, unlike in the numerical studies above, we will not assume /i(x) 
to be spatially uniform. 
At given temperature and chemical potential the expansion then takes the form 

n[M, /2] = Q[0, /io] + -^ / dx OGL(Af (x), (^/i(x)) , (33) 

with 

f^GL(M(x),,5/i(x)) 

= C2,a|M(x)|2 + C2,bS[l{xf + C3,a|M(x)|25/i(x) + C3,bSfi{xf 

+ C4,a|M(x)|^ + C4,b|VM(x)|2 + C4,c|M(x)|2<5A(x)2 + C4,d5/i(x)4 + C4,e(V5/i(x))2 + . . . , (34) 

when expanding to fourth order in M(x), 6jl{x.) and gradients acting on these functions. The 
symmetries of the theory and of the background simplify the expansion: linear terms vanish as we 
are expanding around a homogeneous solution of the gap equations, odd terms in M(x) vanish by 
chiral symmetry and odd numbers of derivatives vanish due to rotational symmetry. 
Taking M(x) to be the small scale of interest, we first aim at an estimate for 5/i(x; M(x)) defined 
through the stationary constraint 



5n 

55jl 



= 0. (35) 

A/(x),5/i(x)=5/i(x;Af(x)) 



Because of the absence of a linear term in M(x) we conclude that 5/i(x;M(x)) ~ 0(|M(x)p). 
More precisely, 

5A(x;M(x)) = -^|M(x)|2 + .... (36) 

Consequently, the expansion of Q.GL{M{ii.)) = CIgl{M{x), (5/i(x; M(x))) to fourth order is given by 

ncLiMix)) = OGL[0,/Uo]+C2,a|M(x)P+ (c4,a-^ ) |M(x)|4 + C4,fe| VM(x)|2 + . . . ,(37) 

V 4c2,6y 

which allows us to determine the Lifshitz and the critical points from the GL coefficients. The 
latter is characterized by vanishing quadratic and quartic mass terms, while at the former the 
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quadratic term and the gradient term are zero: 

LP: = C2,a = C4,fe , 

CP: = C2,a = C4,a " T^ • (38) 

4c2,b 

As outlined in Ref. [44 ] for the NJL model without vector interactions, it is a straightforward 
exercise to work out the explicit form of the GL coefficients. As obvious from Eqs. Q and ([5]), 
the kinetic part of the thermodynamic potential, J^kinotici and, hence, its contributions to the GL 
coefficients depend on Gy only indirectly through JIq and T. The only explicit dependence on Gy 
therefore originates from ilcond and only affects C2,6. For this reason the Lifshitz point as a function 
of JIq and T is independent of Gy , in agreement with our findings in section III CI 
In contrast there is an explicit dependence of the location of the critical point on Gy through C2^b- 
One finds 

C2,6 = -Jt- -nJ~4 + ^], (39) 



4Gy " V^^ 3 

where the first term on the right hand side is due to r^cond) while the second term is due to r^kinetic- 
Since we are expanding around a homogeneous restored solution, the latter is just given by the 
corresponding term in an ideal Fermi gas at temperature T and chemical potential /io- Here we 
have neglected the contributions from the regulator terms, which would arise in our specific model. 
However, these terms are small in the region of interest and therefore do not lead to qualitative 
changes. Hence, l/c2,b vanishes for Gy = and decreases monotonously with increasing positive 
values of Gy. Furthermore c^^a typically decreases when increasing /xq or decreasing T in the 
vicinity of the critical point. Put together, we conclude from Eq. (I38p that the CP moves to 
smaller temperatures upon increasing Gy. Moreover, discarding possible issues related to the 
regularization of the UV-divergent vacuum contribution to the thermodynamic potential^, we find 
C4,a = C4fi, as in the case without vector interaction [44]. For Gy = we therefore recover the 



result of Ref. 4J] that the LP and the CP coincide. 

Since 5/i(x) ~ 0{M{x.)'^) we can also conclude that the thermodynamic potential expanded around 

jlo to order 0(]Af (x)]^) and arbitrary gradients coincides with that of the model for Gy = upon 

replacing // — )• /xq. As a result, the second-order phase transition from any inhomogeneous to the 

chirally restored phase, being triggered by these contributions, is only depending on Gy through 

/io, as it was already obtained in section III CI by applying further truncations. 

On the other hand, the present analysis is not applicable to the transition from the homogeneous 

broken to inhomogeneous phase, where the mass function is not related to a small parameter, when 

we go away from the LP. As we have argued earlier, even the order of the phase transition may 

change in this regime when vector interactions are included. 

Qualitatively the same properties as discussed here for the NJL model with a vector interaction 

also show up in the Gross-Neveu model with a Thirring interaction (GNT model) at leading order 



For a renomializable theory divergent GL coefficients are subject to renormalization; the present discussion is 
consistent for a regularization scheme acting on the energy spectrum as apphed in this work. 
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in the large A^-expansion. Since the model is renornializable and therefore much cleaner and easier 
to handle, we include a discussion of its phase diagram in Appendix 1X1 

III. EFFECTS OF POLYAKOV LOOP DYNAMICS 
A. Inclusion of the Polyakov loop 

In order to mimic features of confinement, in particular to suppress the contribution of free con- 
stituent quarks in the confined phase, and in order to include gluonic contributions to the pressure, 



the NJL model can be coupled to an effective description of the Polyakov loop [13l . [ij]. The 
resulting model is known as the Polyakov-loop extended Nambu-Jona-lasinio (PNJL) model. 
The Polyakov loop is defined by 



L(x) = "Pexp 



i I (irA4(r, x) 
Jo 



(40) 



where A4(r, x) = iAQ{t = —ir,^) is the temporal part of a gauge field A^ = gA^^— at imaginary 
time. In pure Yang-Mills theory, the traced expectation values of L and its hermitean conjugate, 

£=^(T:,L), £-=-l(Tr,Lt), (41) 



can be related to the free energies of a static quark or antiquark, I ~ e~^''i'^ , I ~ e"^''' "^ 50|, |51 |. 
and are therefore order parameters for the confinement-deconfinement transition. 
To include the Polyakov-loop dynamics in the NJL model, the quarks are minimally coupled to a 
background gauge field, d^ — )• D^ = d^ + iAqS^q. Furthermore, a local potential l{{i,£) is added 
to the thermodynamic potential, which is essentially constructed to reproduce ab-initio results of 
pure Yang-Mills theory at finite temperature |14l. |52||. In the most simple approach, the gauge field 
is taken to be a space-time independent mean field A4, and the effect of the covariant derivative 
in the kinetic part amounts to the replacement |13l | 



McrmAE) ^ /thc™ai,PNJL(i?) = Tin (l + e-^^^-^^^ + 3ie-^^-^^/^ + 3^e-2(^-/^)/^; 

+ r In f 1 + e-3(^+'^)/^ + 3 ^e-(^+^)/^ + 3 £ e-2(^+^)/^) . (42) 



However, as can be seen from Eqs. (j40p and ()4ip . a constant mean field ^4 would always result 
in complex conjugate values for i and i. This should be considered as an artifact, since their 
interpretation as being related to the free energies of quarks and antiquarks means that i and i 
should be real and, at finite chemical potential, different from each other. In order to by-pass this 



problem, we therefore follow the viewpoint of Ref. 10| that, after the replacement (j42p . i and i 
should be treated as independent real mean fields, rather than the components of ^4. This is also 
consistent with the fact that the potential U which is added to the thermodynamic potential is 
given in terms of i and i as well. For simplicity we take Fukushima's parametrization [10], 

U{i,e) = -bTlbAe-^/^ei + logll-Qli-^iUf + Aie^+P)]) (43) 
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FIG. 8: Phase diagram of the NJL (sohd hne) and PNJL (dashed Ime) model allowmg for one-dimensional 
spatial modulations of the order parameter. 



with two parameters a and b. Other prescriptions, like the polynomial [IJ] or the logarithmic 52 1 
potential, would lead to very similar results. 

When dealing with inhomogeneous phases, i and £ are naturally expected to be spatially dependent, 
presumably following the density profile in some way. Nevertheless, similar to the treatment of 
fL in the previous section, we will assume spatially independent values of £ and i, even in the 
inhomogeneous phase. This is not only to keep the technical side of the calculation trackable, 
but also the assumptions made in order to derive (j42p and the unknown kinetic contributions to 
Eq. (US]) suggest such a conservative approach as a first step. 
To summarize, we obtain a thermodynamic potential 



Qt 



n^ 



+ ^conA+U{£j), 



(44) 



iPNJL - "kinetic I /^^^^^^i_^/^i^^^^^l_Pl^jL 

with Okinetic and ficond given in Eqs. ()l|)- pT]) . that additionally needs to be extremized in i. and (.. 



B. Numerical results 



Since we consider the NJL model with Gy = ?n- = as our starting point, we shall limit ourselves 
to this case when studyiiig the role of the Polyakov loop. For the Polyakov-loop potential, we adopt 
the parameters of Ref. 10|, a = 664MeV and b = 7.55 • lO^MeV^. The parameter a was fixed 
by the condition that for pure gluo-dynamics the phase transition takes place at T = 270 MeV, 
while b was chosen to have a crossover around T = 200 MeV at ^ = when quarks are included. 
Since we are mainly interested in the qualitative effect of the Polyakov loop, we did not perform a 
refit of b within our regularization scheme. We checked, however, that this parameter choice gives 
reasonable results for the behavior of the order parameters at // = 0. 
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FIG. 9: Polyakov loop expectation value £ in the /i-T plane in the vicinity of the inhomogeneous phase. 



In Fig. [8] we compare the phase diagrams for the NJL model with that of the PNJL model, 
allowing for phases with a one-dimensional solitonic modulation in both cases. Aside from a 
general stretching in the T-direction, which is well known from studies of homogeneous phases 
and easily explained by the replacement (|42p . the Polyakov loop has no effect on the qualitative 
structure of the phase diagram. In particular, the critical point at Gy = still coincides with the 
Lifshitz point. 

In Fig. [9] the value of i is presented via color coding in the region of the phase diagram where the 
inhomogeneous phase is favored. We find that i and i are rather small in the entire inhomogeneous 
phase, reaching their maximum values ^ ~ 0.15 and i ~ 0.2 near the LP. In this context we should 
recall that at vanishing temperature the Polyakov-loop dynamics decouples completely from the 
quark sector due to the way the PNJL model is constructed. As a consequence, ^ = ^ = OatT = 0, 
independent of the density. While it is unclear whether this feature of the model is realistic, it 
means that our assumption of space-independent Polyakov-loop expectation values cannot have a 
large effect. Even if i and i followed the density profile, the results would not be very different, 
because at low temperatures their values are very small anyway, whereas at higher temperatures 
the density differences get washed out. 



IV. DISCUSSION 



In this work we have analyzed the role of the isoscalar- vector interaction and the dynamics of the 
Polyakov loop on recently discussed inhomogeneous ground states in the phase diagram of the NJL 
model. Mainly for technical reasons we thereby limited ourselves to inhomogeneous phases with a 
one-dimensional modulation, explicitly to domain-wall soliton lattices and for comparison to chiral 
spirals. This allowed us to exploit the knowledge obtained for lower dimensional models in our 
study. 
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Our previous study in absence of vector interactions has led us to the finding that the critical point 
in the mean-field phase diagram of the NJL model actually coincides with a Lifshitz point when 
allowing for the possibility of inhomogeneous phases. The first order phase transition is then absent 
in the phase diagram, since an inhomogeneous phase is energetically preferred in its domain. When 
extending our model, we find that a repulsive vector interaction leads to significant qualitative 
effects: In contrast to the critical point when limiting to homogeneous phases, the Lifshitz point 
is not shifted towards smaller temperatures when increasing the strength of the vector interaction, 
but remains at the same temperature and density. Moreover, the domain of inhomogeneous ground 
states in the /U — T phase diagram increases. Since the Lifshitz point therefore no longer coincides 
with the critical point, the critical behavior in its vicinity and also near the phase transition lines to 
the inhomogeneous phase changes. This is underlined by the fact that the number susceptibilities 
remain finite, i.e., there are no singularities. 

Our investigation of this part is complemented by an extensive numerical study, including the 
determination of density profiles in the inhomogeneous phase, the phase diagram when only al- 
lowing for chiral spirals, the evaluation of number susceptibilities and the incorporation of finite 
current quark masses. Furthermore, we performed a generalized Ginzburg-Landau expansion for 
elaborating the dependence of Lifshitz and critical point on the vector interaction and we mapped 
out the phase diagram of the 1 -|- 1-dimensional GN model with a Thirring interaction. Part of this 
comprehensive study is also motivated in order to back up an approximation employed within our 
mean-field calculation, namely to use the spatially averaged number density when evaluating the 
thermodynamic potential. 

Probably less spectacular but nevertheless worth checking is the behavior of our model when the 
quarks are coupled to the Polyakov loop in order to suppress their contribution to the thermody- 
namic potential in chirally broken phases, thus mimicking confinement. In the absence of a vector 
interaction this coupling, at least in the employed approximation, does not lead to a separation 
of Lifshitz and critical point. Consequently the fJ- — T phase diagram is not modified qualitatively 
and, similar as in the case of homogeneous phases, is only stretched in the temperature direction. 
Our analysis has shown that the region where solitonic modulations of the chiral order parameter 
are favored over homogeneous phases increases when considering natural extensions of the two- 
flavor NJL model. Although being much more tractable on the technical side, it is worth noting that 
strictly speaking one-dimensional modulations of the order parameter are washed out by thermal 



fluctuations at finite temperature [53l . |54| | . This statement refers to the fact that the modulation is 
not rigid, but fluctuates locally. The system's spatially modulated nature is however still encoded 
in the behavior of long range correlations. To get past this kind of issue, it would anyway be 
of great interest to extend the present analysis to inhomogeneous phases with higher dimensional 



modulations. This should be possible in a more numerical approach as outlined in Ref. 3]| in 
the context of inhomogeneous color-superconducting phases. The most interesting question here is 
whether the transition from the inhomogeneous to the restored phase can turn first order for more 
complex modulations at low enough temperatures. E.g. for color-superconducting phases, this has 



been suggested in Ref. [271], although the applied expansion did not allow for a conclusive answer. 
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Appendix A: The Gross-Neveu model with a Thirring interaction 

To back up our results related to the vector interaction we briefly want to discuss its effect in the 
Gross-Neveu model with a Thirring interaction (GNT model). Unlike the NJL model, the GNT 
model is renormalizable and therefore does not suffer from regularization artifacts. The relevant 
Lagrangian is given by 

£ = ^ {ij^d^ -m)i; + ^ {^^P) ' - ^ (V^T^V') ' , (Al) 

where ^ is a 2A^-dimensional spinor for N species in 1 + 1 dimensions and 7^* can be chosen as 
7^ = (T^, 7^ = —ia"^, 0"* being the usual Pauli matrices. 

In the large A^-limit the mean-field approximation becomes exact. ^ Furthermore, since the model 
is renormalizable, all divergencies can be absorbed into a finite number of contact terms. The 
thermodynamic potential is given by 

^GNt/^ = ^kinetic /A^ + ^cond/^ , 

^kinetic /A^ = --^TrD,LLogf-(iw„ + F)j , 

.,„,/^,.l/..(M£l_:!i)!_M£^_il)!), (A2) 

where we introduced a similar notation as in the case of the NJL model: L is the periodicity, the 
spatial coordinate is labelled z, M{z) = m — Gs{'>l'ip)/N, jl{z) = /U — Gvi'4'^^tp)/^ aiid 

F = -i7%ia, + 7°M(z) - /i(z) . (A3) 

For this setup we first investigate the phase diagram of homogeneous phases, in particular the role 
of Gy, and subsequently also inhomogeneous phases by means of a Ginzburg-Landau expansion. 
For homogeneous phases the eigenvalue spectrum of F can be labelled by the momenta and it is 
straightforward to obtain 



no^T/N = -2/ £e,+ '—--^-'J^-£L_2T I £T.'n(l + e-^, (A4) 



s=± 



— ^'vacuum , ,^ 



'-vacuum 



Note that this hmit does not commute with the thermodynamic hmit. 
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where Ep = yp^ + M^ . For r^vacuum we introduced a sharp momentum cutoff A as the integral is 
divergent. Note, however, that for a renormahzable theory - in contrast to the NJL model - the 
regularization scheme is only introduced intermediately and does not affect the final results. We 
also note that the argument of the vacuum contribution should in principle contain the eigenvalues 
±Ep — JJ, + jjL,'^ but for a sharp momentum cutoff this does not affect the renormalization procedure 
at large momenta. Concerning the ultraviolet behavior, we observe that the quadratic divergency 
can be absorbed into an order parameter independent constant, a linear and a quadratic term in 
M, corresponding to a renormalization of the overall pressure at M = 0, the fermion mass m and 
the coupling Gs- Requiring the thermodynamic potential to vanish at M = {ip^^ip) = and to 
have a minimum at Mq, we then get 

vacuum/ ^^ ^ ^ Mq J 2Mq 2Gv 

where c is proportional to m at finite A and not relevant here as we shall consider the chiral limit 
c = in the following. Gy remains finite and can be assumed to have its renormalized value. In 
order to discuss the phase diagram all that is left is to extremize the potential in M and fl. For 
simplicity we will focus the discussion on phase transitions and the location of the critical point. 
Addressing phase transitions between homogeneous phases, we first expand the thermodynamic 



potential in the constituent mass around M = 0. Similar to the GN model 55| we find 



^gnt/N = ^a2„(/i)M2" (A6) 



n>0 



with 



ao(/i) 



ttT' fi' (/i - nY 



6 27r 2Gv 



(-1)" 

where we used the di- and polygamma functions and z = 2 + 2^ • Note that only 02 is affected by 
the renormalization procedure. From the additional constraint dVLQ^-^ / djl = we can then infer 
the value /2o = ^^q of the renormalized chemical potential at M = and furthermore expand in 
5jl = jl — JIq. For the purpose of our discussion, we can then limit ourselves to 

^gnt/A^ = ^unbroken + a2fiM'^ + a^fiM* + a2,iM'^5fL + Qo,2<5/i^ + 0{M\ 5jj.\ M^6jl^, M^5fL) , 

(A8) 



This is related to the fact that, Uke the constituent mass M, the density {ijj'y'^ip) oc p — /i is undetermined before 
the gap equations are solved. 
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where an,o = an|/i=/io and 

^ ^unbroken = CKO.O 



7rr2 fi 



2 



6 2(7r + Gy)' 
1 



02,1 = -^;f327^ImV'(l,z)|^=^o' 



TT + Gy 

We conclude that (5/x = — 2a'^ -^^ + O(M^) and arrive at an expansion of the thermodynamic 
potential only depending on M, 

Qgnt/N = f^unbrokcn + "2,0^2 + /34M4 + 0(M6) , (AlO) 

where 

/34 = 04,0 - 7^^ • (All) 

The critical point is then given by 02,0 = /34 = 0. As a simple exercise we can evaluate this 
condition at zero temperature, for which one finds ln(2/io/Mo) = — 1 + -Gy = 0. At Gy = 7r/3 
we have therefore a critical point at T = and /2 = 2Mq/3. However, as illustrated in Fig. ijlOp . it 
turns out that this is a "new" CP, which moves upwards in temperature upon increasing Gy- On 
the other hand, the "old" CP, existing already at Gy = 0, moves downwards, similar to the NJL 
model. Eventually, at Gy ~ 1.175, both points merge and no first-order phase transition is left at 
higher values of Gy . 

Turning to inhomogeneous phases, we can perform a generalized Ginzburg-Landau expansion. As 
in the NJL-model case, discussed in Sec. IIIDl this corresponds to an expansion of ()A2p in M(z) 
and 5fl{z) = jl{z) — (1q, combined with a derivative expansion in order to obtain a local functional. 



The technology has been worked out in great detail for the GN model 3a| and equally works 
for the GNT model. Since a M{z)6p,{z)-tej:m is forbidden by Z2 symmetry, we have in general 
6fi{z) ~ M(z)2. Furthermore, by treating derivatives to be of order 0{M{z)) we get to fourth 
order 

^GNT/iV = f^unbrokon + -^ dz(^a2,oM{zf + a4,o{M{z)^ + M' [zf] + a2,iM^5fl{z) 

+ao,26fx{zf) + . . . (A12) 



and observe that the prefactors of the M{z)'^ and M'{z)'^ terms are those obtained for the GN 
model upon replacing // — )• /lo • From this we conclude that the location of the Lifshitz point in the 
fl — T diagram is the same as in the GN model. Consequently, since the density there is directly 
given through [Iq, Gy does not affect the Lifshitz point in the n — T diagram, while in the fi — T 
diagram the LP is shifted by Gy{'ij}'y^ip)/N in the //-direction. 

Another property that is inherited from the GN model is the second-order transition line from the 
chiral crystalline to the restored phase. Since in this case the magnitude of the order parameter 
vanishes continuously while the modulation stays finite, we have to consider the GL expansion to 
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FIG. 10: Phase diagram of the massless GNT model: The dotted line corresponds to q;2.o — 0, the dots 
indicate the critical points when limiting to homogeneous phases, i.e., in addition /34 = 0. At Gy = 7r/3 
an additional CP emerges on the /i-axis, which merges with the other CP at Gv,max — 1.175. The solid 
line shows the transition from the chiral crystalline to the restored phase, ending at the Lifshitz point. The 
latter agrees with the CP at Gy — 0. 



order 0{M{z)'^) and in principle to all orders of gradients. However, since Sfl{z) ~ M(z)^, 6fl does 
not enter in this analysis and we have to have the same result as in the GN model. For the latter 
this transition line has been determined in Refs. 3g, |38|] and is also depicted in Fig. (|10p . 
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